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I. INTRODUCTION 

Real-time methods have proven their utility in calcu- 
lating optical properties of finite systems mainly within 
time-dependent density functional theory (TDDFT)iii^ 
On the other hand extended systems have been 
mostly studied by using many-body perturbation theory 
(MBPT) within the linear response regime^ The differ- 
ent treatment of correlation and nonlinear effects mark 
the range of applicability of the two approaches. The 
real-time TDDFT makes possible to investigate nonlinear 
effects like second harmonic generatiorki or hyperpolar- 
izabilities of molecular systemsi^ However the standard 
approaches used to approximate the exchange-correlation 
functional of TDDFT treat correlation effects only on a 
mean-field level. As a consequence, while finite systems — 
such as molecules — are well described, in the case of 
extended systems — such as periodic crystals and nano- 
structures — the real-time TDDFT docs not capture the 
essential features of the optical absorption^ even qualita- 
tively. 

On the contrary MBPT allows to include correlation 
effects using controllable and systematic approximations 
for the self-energy E, that is a one-particle operator non- 
local in space and time. E can be evaluated within dif- 
ferent approximations, among which one of the most suc- 
cessful is the so-called GW approximatioui^ Since its first 
application to semiconductors^i the GW self-energy has 
been shown to correctly reproduce quasi-particle energies 
and lifetimes for a wide range of materials Furthermore, 
by using the static limit of the GW self-energy as scat- 
tering potential of the Bethe-Salpeter equation (BSE)^ 
it is possible to calculate response functions including 
electron-hole interaction effects. 

In recent years, the MBPT approach has been merged 
with density-functional theory (DFT) by using the Kohn- 



Sham Hamiltonian as zeroth-order term in the perturba- 
tive expansion of the interacting Green's functions. This 
approach is parameter free and completely ah-initio^ and 
in this paper will be addressed as ah-initio-MSPT {Ai- 
MBPT) to mark the difference with the conventional 
MBPT. However the 4j-MBPT is a very cumbersome 
technique that, based on a perturbative concept, in- 
creases its level of complexity with the order of the ex- 
pansion. As an example, this makes the extension of this 
approach beyond the linear response regime quite com- 
plex, though there have been recently some applications 
of the ^i-MBPT in nonlinear optics<^""— 

Another stringent restriction of the ^i-MBPT is that it 
cannot be applied when non-equilibrium phenomena take 
place: for example it cannot be applied to study the light 
emission after an ultra-fast laser pulse excitation. A gen- 
eralization of MBPT to non-equilibrium situations has 
been proposed by Kadanoff and BaymJ^ In their semi- 
nal works the authors derived a set of equations for the 
real-time Green's functions, the Kadanoff-Baym equa- 
tions (KBE's), that provide the basic tools of the non- 
equilibrium Green's Function theory and allow essential 
advances in non equilibrium statistical mechanics 

Both the standard MBPT and non-equilibrium Green's 
Function theory are based on the Green's function con- 
cept. This function describes the time propagation of a 
single particle excitation under the action of an exter- 
nal perturbation. In the equilibrium MBPT, due to the 
time translation invariance, the relevant variable used to 
calculate the Green's functions is the frequency lj. In- 
stead, out of equilibrium, in all non steady-state situa- 
tions, the time variables acquire a special role and much 
more attention is devoted to the their propagation prop- 
erties. The time propagation avoids the explosive de- 
pendence, beyond the linear response, of the MBPT on 
high order Green's functions. Moreover the KBE are 
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non-perturbative in the external field therefore weak and 
strong fields can be treated on the same footing. 

One of the first attempts to apply the KBE's for in- 
vestigating optical properties of semiconductors was pre- 
sented in the seminal paper of Schmitt-Rink and co- 
workers!^ Later the KBE's were applied to study quan- 
tum wells laser excited semiconductors)^ and lumi- 
nescencei^. However, only recently it was possible to 
simulate the Kadanoff-Baym dynamics in real-timeJ^"— 

In this work we combine a simplified version of the 
KBE's with DFT in such a way to obtain a parameter- 
free theory that is able to reproduce and predict ultra- 
fast and nonlinear phenomena (Sec. This approach, 
that we will address as time-dependent BSE, reduces to 
the standard BSE for weak perturbations (Sec. Ill C|) but, 
at the same time, naturally describes optical excitations 
beyond the linear regime. After discussing some relevant 
aspects of the practical implementation of our approach 
(Sec. imp , we exemplify how it works in practice by cal- 
culating the optical absorption spectra of /i-BN and the 
time dependent change in its electronic population due 
to the perturbation by means of an ultra-fast and ultra- 
strong laser pulse (Sec. lIVp . 

II. THE TIME-DEPENDENT 
BETHE-SALPETER EQUATION 

We derive here a novel approach to solve the time evo- 
lution of an electronic system with Hamiltonian coupled 
with an external field, 

H ^h + Hrab + U, (1) 

where U represents the electron-light interactions (see 
Sec. IIII Al for its specific form). As usually done in 
MBPT, H is partitioned in an (effective) one-particle 
Hamiltonian h and a part containing the many-particle 
effects H„ib- 

In our derivation, we take as starting point the KBE's 
that we briefly introduce in Sec. Ill Al fsee e.g. Refs.[20lfor 
a systematic treatment). Then, in Sec. Ill Bi we proceed in 
analogy with the equilibrium j4i-MBPT: first, we define 
h as the Hamiltonian of the Kohn-Sham system, second 
we introduce the same approximations for the self-energy 
operator. As a result we obtain the analogous of the 
successful Giy-|-BSE approach for the non-equilibrium 
case. Indeed in Sec. Ill CI we show that our approach, 
the time-dependent BSE, reduce to the GM^-I-BSE in the 
linear regime. 

A. The Kadanoff-Baym equations 

Within the KBE's, the time evolution of an electronic 
system coupled with an external field is described by 
the equation of motion for the non-equilibrium Green's 
functionsiii^Si^i, G (r, i; r'i'). To keep the formulation 



as simple as possible and, being interested only in long 
wavelength perturbations, we expand the generic G in 
the eigenstates {fn,\^] of the fi Hamiltonian for a fixed 
momentum point k: 

[Gk (ii, i2)]„j„2 = G„i„2,k(ii, ^2) = 

j ^;^k(i'i)G(ri,ii;r2,i2)¥'n2k(r2)dVidV2. (2) 

As the external field does not break the spatial invariance 
of the system k is conserved. 

Within a second-quantization formulation of the many- 
body problem, the equation of motion for the Green's 
function described by Eq. ^ are obtained from those 
for the creation and destruction operators. However the 
resulting equations of motion for Gk are not closed: they 
depend on the equations of the two-particle Green's func- 
tion that in turns depends on the three-particle Green's 
function and so on. In order to truncate this hierar- 
chy of equations, one introduces the self-energy opera- 
tor 5]k(ii,t2), a non-local and frequency dependent one- 
particle operator that holds information of all higher or- 
der Green's functions. A further complication arises with 
respect to the equilibrium case because of the lack of 
time-traslation invariance in non-equilibrium phenomena 
that implies that Sk(^i, ^2) and Gk(ti, ^2) depend explic- 
itly on both ti,t2- Then, one can define an advanced 
(G^), a retarded SJ^ (G^), a greater and a lesser S^, 
(G^,G^) self-energy operators (Green's functions) de- 
pending on the ordering of ii,i2 on the time axis. Fi- 
nally, the following equation for the is obtained (see 
e.g. Ch. 2 of Ref. [20 for more details): 

d 

«^^G'< „^k(^l' ^2) =5{tl- t2)Sn,n, 

(^1,^2) +^ C/ni«3k(il)G<3„2k(^l>*2) 

"3 

+ E / d^3(S^Hn3k(^l:i3)G<„^kfe,i2) 
na 

+ ^rHn3k(^l'*3)G^3„2k(^3,^2))- (3) 

This equation, together with the adjoint one for ih-^C^ , 
describes the evolution of the lesser Green's function G^ 
that gives access to the electron distribution {G^{t,t)) 
and to the average of any one-particle operator such as for 
example the electron density [Eq. (fTO)) ]. the polarization 
[Eq. p2p ] and the current. However, in general E'',E< 
and the G^ depend on G^ , so that in addition to Eq. ([3]) 
the corresponding equation for the G^ has to be solved. 

Then, in principle, to determine the non-equilibrium 
Green's function in presence of an external perturba- 
tion one needs to solve the system of coupled equations 
for G^,G^, known as KBE's. Indeed, this system has 
been implemented within several approximations for the 
self-energy in model systems,— in the homogeneous 
electron gas,— and in atoms^. The possibility of a di- 
rect propagation in time of the KBE's provided, in these 
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systems, valuable insights on the real-time dynamics of 
the electronic excitations, as their lifetime and transient 
effects.—"— Nevertheless, the enormous computational 
load connected to the large number of degrees of free- 
dom de facto prevented the application of this method 
to crystalline solids, large molecules and nano-structures. 
In the next subsection we show a simplified approach — 
grounded on the DFT — that while capturing most of the 
physical effects we are interest in, makes calculation of 
"real- world" systems feasible. 



B. The Kohn-Sham Hamiltonian and an 
approximation for the self-energy 

In analogy to j4i-MBPT for the equilibrium case, we 
choose as h in Eq. ([T|) the Kohn-Sham Hamiltonian,— 



2m ^ 



Vi + Vei + V"[~p]+V^^[p\, (4) 



where Vei is the electron-ion interaction, is the 

Hartree potential and F^'^ the exchange-correlation po- 
tential. Within DFT, the Kohn-Sham Hamiltonian cor- 
responds to the independent particle system that repro- 
duces the ground-state electronic density p of the full 
interacting system (ft, -I- -ffmb), that is 



E 

nk 



/nk|^(r)| = 



(5) 



where /„k is the Kohn-Sham Fermi distribution. 

Equation ([3]) can be greatly simplified by choosing a 
static retarded approximation for the self-energy. 



S<(ii,i2) =0 



(6a) 
(6b) 



where the usual choice is S^ohscx^ ^^^^ so-called Coulomb- 
hole plus screened-exchange self-energy (COHSEX). In 
Eq. |6a] we subtracted the correlation effects already ac- 
counted by Kohn-Sham Hamiltonian The COHSEX 
is composed of two parts: 



I]"°'^(r,r',t) = iM^(r,r';G<)G<(r,r',t), 



(7) 



S™'Xr,r',i) 



1 



-M^(r,r';G<)-(5(r-r'), (8) 



where VF(r,r';G^) is the Coulomb interaction in the 
random-phase approximation (RPA). These two terms 
are obtained as a static limit of the GW self-energy (see 
Ch.4 of Ref. HO and Refs. [H and [H . 

With the approximation in Eqs. ([5a|) - (|6b)) . Eq. (jS]) docs 
not depend anymore on and it is diagonal in time: 

*ft|G< ,„^,,(i) = 

[hk + Uk(t)+V^[p]-V^[p] 
^^5,cohsc.(^)_Vr[/5]),G<(t)]„^^„^. (9) 



where p is the density obtained from the G^ as 

nin2k 

Equation ([9]) is conservingii and satisfies the sum rules 
for the response functions because both the one- and two- 
particle self-energies are obtained from the same 
and the system of equations is solved self-consistently<^ 

However, despite the full real-time COHSEX dynam- 
ics [Eq. |9] is an appealing option considerably simplify- 
ing the dynamics with respect to the KBE's, it neglects 
the dynamical dependence of the self-energy operator. 
This, in practice, induces a consistent renormalization of 
the quasiparticle chargo^ in addition to an opposite en- 
hancement of the optical properties^. In the COHSEX 
approximation both effects are neglected. At the level of 
response properties for most of the extended systems dy- 
namical effects are either negligible or very small (while 
recently it has been shown their importance for finite sys- 
tems, see Refs. [13 and HI) and, for practical purposes, it 
has been shown that they partially cancel with the quasi- 
particle renormalization factors.— 

Therefore we modify Eq. IH] in order to include only 
the effect of the dynamical self-energy on the renormal- 
ization of the quasi-particle energies, that is the most 
important effect. Also in this case, our idea is to pro- 
ceed in strict analogy with ^?-MBPT and to derive a 
real-time equation that reproduces the fruitful combina- 
tion of the Go Wo approximation — for the one-particle 
Green's function — and of the BSE with a static self- 
energy — for the two-particle Green's function. Indeed 
the GoWo-f BSE is the state-of-the-art approach to study 
optical properties within the ^i-MBPTi^ To this purpose 
Eq. ^ is modified as: 



[hk + Ahk + Uk + [p] - [p] 



^cohscx \/2< 



[G<],G< it) 



(11) 



Ah is a scissor operator— that applies the Go Wo correc- 
tion to the Kohn-Sham eigenvalues, e^^^, 



[Ah, 



,GoWo 



^riik) ^ni,n2 7 (12) 



and G^„, is the solution of Eq. ([TT|) for the unperturbed 
system ([/ = 0) 



G 



< 

nn'k 



^^y'nk'^nr: 



(13) 



where we assume that the Kohn-Sham Fermi distribution 
is not changed by the scissor operator. Note further that 
in Eq. [TTl V^'^ [p] cancels out because it is independent of 
G<it). ^ ^ 

Equation is the key result of this work. It is equiv- 
alent to assume that the quasi-particle corrections mod- 
ify only the single particle eigenvalues leaving unchanged 



4 




integration of time-dependent BSE 
by 2nd order Runge-Kutta: 

G<{t) 




; Post-processing 

optical properties 



e.g. X,X^^\X^^^ 



FIG. 1. Schematic flow of a time-Dependent BSE simulation. See 
Sec. IIII Xl for details. 



Wc start by expanding x(r) in terms of the Kohn-Sham 
orbitals: 



Z,m,k 



X V3,,k(r)v3*k+q(r)v?r,k'(r')¥'rn,k'+q(r'), (16) 

where q is the momentum, and we define the matrix el- 
ements of x' a-s, 



d^rdVV* k(r)^:„,k'+q(r')^j.k+q(r)¥.i,k' (r')- (17) 



Since we are interested only in the optical response, in 
what follows we restrict ourselves to the case q = and 
drop the q dependence of y'' (for the extension to finite 
momentum transfer see Ref. llSf ) . Inserting the expansion 
for X [Eq. (USD], p [Eq. ^] and U (C/™,„k = {mk\U\nk}) 
in Eq. ([15]) we obtain the following relation linking the 
matrix elements of x"" to the matrix elements of G"^ : 



the Kohn-Sham wave functions. Within ^iMBPT this 
approximation is very successful for a wide range of ma- 
terials characterized by weak correlations (see e.g. Refs.H 
and@). 



Im.p 



SUira.pit') 



(18) 



(7=0 



Then, we can obtain the equation of motion for the ma- 
trix elements of x"^ by taking the functional derivative of 
Eq. (fTT|) with respect to J7;,,n,k {t), 



C. The linear response limit 

When an external perturbation U (t) is switched on in 
Eq. (fTTj) . it induces a variation of the Green's function, 
AG^(t) = G^(t) — G^. In turns, this variation induces 
a change in the self-energy and in the Hartree potential. 
In the case of a strong applied laser field these changes 
depend on all possible orders in the external field. How- 
ever for weak fields the linear term is dominant. In this 
regime it is possible to show analytically that Eq. (fTTj) 
reduces to the GqWo+BSE approach^i^. Proceeding sim- 
ilarly to Ref. [13 we consider the retarded density-density 
correlation function: 



X''{r,t;r',t')^-^[{p{r,t)p{r',t')) 

~{p{r,t)){p{r',t'))]e{t^t'). (14) 



d 

^'^ lm,p 

— ^— -[hk + Ahk + Uk(t) + [p{t)] - [p] 



Sk[G<(t)]-Sk[G<],G<(t)] 



(19) 



Making use of the definitions in Eqs. ([T2|) and (|T3|) . to- 
gether with Eq. (fTS)) , it can be verified that the functional 
derivative of the one-electron Hamiltonian and of the ex- 
ternal field give the contribution 



-[hk + Ahk + Uk,G<(t)] 



(7=0 



(e 



Go Wo 



)x',i,k (t - t') + i{fiV - /jk)'5j7'5.™5kp. 



(20) 



X' describes the linear response of the system to a weak 
perturbation, represented in Eq. ([1]) by [/, 



X'{^,t-v\t') 



(Spirt)) 



6U{r't') 



(15) 



Note that, since the perturbation is weak, the Hamil- 
tonian of the system is invariant with respect to time 
translation and thus x^ depends only ont — t'. The term 
in Eq. (fTO)) containing the Hartree potential, that is not 
directly depending on the external perturbation, is ex- 
panded with respect to Ui^mM (t) by using the functional 
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derivative chain rule and the definition in Eq. (|18p as: 



Ti,n ,p 



n'n,p 



xxL.',p(i''^"MC/;m,k'(t"), (21) 



A similar equation can be obtained for Ti1'y^''^[G'^{t)]. 
Equation (PT|) for Hartree potential and its analogous for 
the self-energy can be explicited by using 



q=0 
rin,k 



ji,(k-q) "'mk,i(k-q): 
y,q nk,i(k-q) 



(22) 



(23) 



where the matrix elements of u^^" and W are labeled 
accordingly to Eq. ((TT)) . In Eq. (l22t w ''^^ is the long range 
part of the bare Coulomb potential, responsible for the 
local field effects in the BSE. Then by inserting Eq. ((^ 
in Eq. ([?!]) the functional derivative for the Hartree term 
is 



7Y[Vf[p(0]-V^[p],G<(0]^, 



u=o 



(2»2) (/,k - /,k) E ^J^^x'stM (t ~ t'). (24) 



Similarly, an analogous equation is obtained for the 
self-energy (see also Appendix O, 



u=o 



(/.k - ./,k) E W^.k,s(k-q)X'.t,(k-q)(i - i'), (25) 

st,q lk,t(k-q) im,p 

where we neglected the part containing the functional 
derivative of the screened interaction with respect to the 
external perturbation. This is a basic assumption of the 
standard BSE that is introduced in order to neglect high 
order vertex corrections.— 

Finally, we insert Eqs. ([24]) and ([25]) in Eq. (fT9| . 
and by Fourier transforming with respect to (t — t') we 
obtain 



hio — 



V Jk 



Wo _ ^GoWo 



/?np 



-iE{^ik,s(k-q) 
si,q ik,t(k-q) 



2v 



stM lm,p 



(26) 



III. OPTICAL PROPERTIES FROM A 
TIME-DEPENDENT APPROACH 

A. Practical solution of the time-dependent BSE 

To solve Eq. (fTTj) for a given electronic system [Eq. ([T])] , 
we start from h, with its eigenvalues and eigenstates de- 
termined from a previous DFT calculations, and from the 
corrections Ahk, determined e.g. from a previous GqWq 
calculation. Then, we switch on the external perturba- 
tion U and integrate the equations of motion using the 
same scheme as in Ref. [iB for the diagonal part of the 
G^, that is equivalent to a second order Runge-Kutta. 
Specifically, in Eq. Q we choose to treat the interac- 
tion with the external electric field E within the direct 
coupling — or length gauge. 



(27) 



Other choices are possible and indeed in the literature 
the electron-light interaction is often described within the 
minimal coupling — or velocity gauge (p • A) , with A the 
vector potential. As it has been pointed out in Ref. [2^, 
the length and velocity gauges lead to the same results 
only if a gauge transformation is correctly applied. How- 
ever, in this respect the velocity gauge presents two main 
drawbacks. First, the wave functions and the boundary 
conditions have to be transformed by a time-dependent 
gauge factor T(r, t) ~ exp{iA(t) • f } and accordingly, 
in the Green's function formalism also the self-energy 
and the dephasing term have to be transformed. Second, 
within perturbation theory the velocity gauge induces di- 
vergent terms in the response function that in principle 
cancel each other, but that in practice lead to artificial 
divergences in the optical respons o^'^" due to numerical 
precision and incomplete basis sets. 

The interaction Hamiltonian U is evaluated in terms 
of unperturbed Kohn-Sham eigcnfunctions as 

(mk|J7|nk) = -E(i)(mk|r|nk) = -E(t) r„„^k, (28) 

where the dipole matrix elements rmn,k , for m ^ n are 
calculated by using the commutation relation i[H,r] = 
p + i[Vni, r] where Vni is the non-local part of the Hamil- 
tonian operatori^ 

Since we are interested in calculating the dielectric 
properties at zero momentum, we choose to work with 
an homogeneous electric field E(i), with no space depen- 
dence except its direction^i generated by a vector poten- 
tial A(i) constant in space. 



1 dA{t) 
c dt 



(29) 



formally equivalent to the standard BSE. 



Also in this case other choices consistent with the peri- 
odic boundary conditions would be possible, as for ex- 
ample an external potential with the cell periodicity;^ 
or an electric field with a finite momcnturoi^ q such that 
q = k-k'. 
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Time (fs) Energy (eV) 




Time (fs) Energy (eV) 

FIG. 2. /i-BN: Comparison between the real-time approach and the standard RPA and BSE approaches based on the equihbrium 
MBPT. (a),(c): polarization F{t) generated by an electric field E(t) = EoS{t) within the TD-HARTREE [(a)] and TD-BSE [(c)] 
approximations. (b),(d): the corresponding absorption spectra (red circles) are compared with the RPA [(b)] and with the BSE [(d)] 
results (black line). The experimental absorption spectrum (grey shadow) is also shown as reference. 



Instead, the particular form of E(t) as function of time 
is not specified a priori, but given as input parameter of 
the simulation. Indeed, the possibility of providing the 
form of the external field as an input is one of the key 
strengths of the real-time approach, potentially allow- 
ing to use the same implementation to simulate a broad 
range of phenomena and of experimental techniques. For 
example, as described in Sec. IIVI in order to calculate 
the linear optical susceptibility spectrum xi^) '^^ will 
use a delta function E(i) = Eo(5(t — to) (obtained from 
Eq.dini) with A(t) = Ao9(t - to), where to is the time 
at which the external field is switched on). This elec- 
tric field probes the system at all frequencies with the 
same intensity. Also, in the other example described 
in Sec. IIVI we can use a quasi-monochromatic source 
E(i) = Eg sincjoi exp (-(5^(t — to)^/2) to selectively ex- 
cite the system at a given frequency wq. Furthermore, 
two or more electric fields can be used to simulate e.g. 
pump-probe, sum-of-frequency or wave-mixing experi- 
ments. 

The macroscopic quantity that is calculated at the end 
of the real-time simulation is the induced polarization 
P{t), related to the electric displacement D(r,t) and the 



electric field E(r, t) by the so called material equation: 

D(r,t) = eoE(r,t)+P(r,i), (30) 

that stems directly from the IVIaxwell equations. P{t) is 
obtained from G< [Eq. (HI])] by, 

P(^) = 5Z ''mn,kG'^m,k(0' (31) 

n,m,k: 

and from this quantity we can obtain the optical proper- 
ties of the system under study. 

For instance, within linear response, the electric dis- 
placement D(r,i) is directly proportional, in frequency 
space, to the electric field as D(cj) = e(cj)eoE(a;). There- 
fore the polarization can be expressed as: 

F{uj) = eo(e(c^) - /)E(c^) (32) 

and accordingly the optical susceptibility that describes 
the linear response of the system to a perturbation is 
x(a;) = e{uj)—I. Then, the optical susceptibility x{^) can 
be calculated by Fourier transforming the macroscopic 
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polarization P(w) (or alternatively the current density 
by means of Eq. ([5^ as: 



eoE(a;) 



(33) 



Note that by choosing a delta-like E(t), the Fourier trans- 
form of P(t) provides directly the full spectrum of the 
optical susceptibility x{(jj)- Beyond the linear regime, 



(2) v(3) 



can be ob- 



higher order response functions, x S X 
tained (to calculate e.g. the second- or third-harmonic 
generation) by using a (quasi)monochromatic field as in 
e.g. Ref. IJ; non-perturbative phenomena, such as high- 
harmonic generation, can be analyzed instead from the 
power spectrum (|P(w)|^). 

To summarize, the schematic flow of a time-dependent 
BSE simulation is shown in Fig. [1] as has been imple- 
mented in the development version of thcYAMBO code.— 
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FIG. 3. A,-BN: Percentage of valence electrons pumped to the con- 
duction bands (A^c) by a quasi- monochromatic pulse as a function 
of the fluence. The pulse is centered either at 5.65 eV (blue boxes) 
or 8.1 eV (green circles) calculated within the TD-BSE (black line) 
and the td-HARTREE (red dashed line) approximations. In either 
case, each point corresponds to a separate simulation and the lines 
are drawn to help guide the eye. The inset shows the absorption 
spectra within the TD-BSE (blaclc line) and TD-HARTREE (red 
dashed line) with the arrows pointing at the pump frequencies. 



B. Dissipative effects 

In an excited electronic system dissipative effects are 
present due to inelastic electron scattering and (quasi)- 
elastic scattering processes with other degrees of freedom, 
such as defects or phonons. Both effects contribute to 



the relaxation and decay of excited electronic population 
as well as of the decay of phase coherence, that is to a 
finite dephasing rate. Our approach, Eq. does not 

account for dissipative effects: on the one hand the COH- 
SEX self-energy is real, so that the excitations lifetimes 
are infinite, on the other hand the electronic systems is 
perfectly isolated [Eq. ([T])], so that there is no dephasing 
due to interaction with other degrees of freedom. 

In practical calculations then we introduce a phe- 
nomenological damping to simulate dissipative effects. 
We implemented two different approaches. An a pos- 
teriori treatment, where at the end of the simulation (in 
the post-processing block of Fig. [T]) the polarization (and 
the electric field) are multiplied by a decaying exponen- 
tial function, e"*/"^, where r is an empirical parameter. 
This parameter, that is compatible with the simulation 
length, effectively simulates the dephasing and introduces 
a Lorentzian broadening in the resulting absorption spec- 
trum. This is in the same spirit of the Lorentzian broad- 
ening introduced in the linear response treatment to sim- 
ulate the experimental optical spectra, and has the ad- 
vantage of producing spectra with different broadening 
from the same real-time simulation. Nevertheless this 
approach is limited to the linear response case. 

In order to treat dissipative effects beyond the linear 
regime, an imaginary term is added to the self-energy 

in the form of an additional term Fk 

pearing on the r.h.s. of Eq. (fTTj) . with: 



ap- 



^ ^ )iik 



712^ 



riik""l,"2 



(34) 



where F^"^ and F^''^ are respectively the lifetime of the 
perturbed electronic population and the dephasing rate, 
and are given as input parameters of the simulation. 



IV. EXAMPLES 

To illustrate and validate the time-dependent BSE ap- 
proach and our numerical implementation, we present 
two examples on /i-BN. This is a wide gap insulator 
whose optical properties are strongly renormalized by ex- 
citonic effects and for which all the parameters necessary 
in DFT, GqWq and response calculations?^ are known 
from previous studies »^ 

In these examples we used Eq. (jlip . with and with- 
out including the self-energy term. We refer to the for- 
mer approximation as TD-BSE, and to the latter as TD- 
HARTREE. Within equilibrium MBPT these two ap- 
proximations would correspond to the BSE and RPA, 
and in fact they reduce to BSE and RPA within the lin- 
ear response limit (Sec. Ill C|) . 

In the first example (Fig. [5]), we simulated /i-BN inter- 
acting with a weak delta-like laser fieldi^ As explained 
in Sec. IIII Al a delta-like laser field probes all frequencies 
of the system and the Fourier transform of the macro- 
scopic polarizability provides directly the susceptibility 
[Eq. and thus the dielectric constant [Eq. ([5^ ]. 




FIG. 4. /i-BN; Percentage as function of time of valence electrons pumped to the conduction bands {Nc) by a quasi-monochromatic pulse 
with intensity 10^ kW/cm'^ centered either at 5.65 eV (blue boxes) or 8.1 eV (green circles) and calculated either within the TD-BSE [(a)] 
or the TD-HARTREE [(b)] approximations. The brown shadow represents the fluence as function of time. The inset shows the absorption 
spectra within the TD-BSE (black line) and TD-HARTREE (red dashed line) with the arrows pointing at the pump frequencies. 



Since we use a weak field, we expect negligible nonlinear 
effects. Then accordingly with Sec. Ill C[ the results from 
TD-BSE and TD-HARTREE can be directly compared 
with the BSE and RPA within the standard ^z-MBPT 
approach. Indeed, in Figs.[2](b),[2](d) the imaginary part 
of the dielectric constant (optical absorption) obtained by 
Fourier transform of the polarization in Figs.[2](a),[2](c) is 
indistinguishable from that obtained within equilibrium 
j4i-MBPT, validating our numerical implementation. 

In the second example (Figs. [3]|3]) we exploit the po- 
tentiality of the TD-BSE approach by going beyond the 
linear regime and using a strong quasi-monochromatic 
laser field (see Sec. IIII A]) . This field excites the system 
selectively at one given frequency, moreover it is strong 
enough to induce changes in the electronic population of 
the system. To track these changes, during the dynamics 
we followed the evolution of Nc{%), that is the percent- 
age of valence electrons that are pumped by the electric 
field in the conduction bands (in our simulation we have 
16 valence electrons in the /i-BN unit cell, since core elec- 
trons arc accounted using pseudopotentials) . The total 
number of valence electrons in the system is given by the 
trace of while Nc{t) = —i J2c]i ^cck^^) where c labels 
the empty states in the unperturbed system. 

We performed the simulations^ for different intensi- 
ties of the field (from lO^kW/cm to lO^kW/cm) and 
for two values of the field frequency, 5.65 eV and 8.1 cV, 



that depending on the level of the theory, are either at 
resonance or off-resonance with the system characteristic 
frequencies. More precisely, within TD-BSE 5.65 eV cor- 
responds to the strong excitonic feature in the absorption 
spectrum, while at 8.1 eV the absorption is negligible; 
conversely within RPA at 5.65 eV the absorption is neg- 
ligible, while 8.1 eV corresponds to the strongest feature 
in the spectrum (see inset of Figs.[3]|31). The results of the 
various simulations are summarized in Fig. |3] that shows 
Nc{%) as a function of the fluence — the pulse energy per 
unit area. For a comparison the ablation threshold of h- 
BN has been determined as 78mJ/cm^ in the femtosec- 
ond laser operational mode<^ 

Finally, Figs. HKa) and HKb) show the evolution of 
Nc{%) during the simulation for a field intensity of 
lO^kW/cm: one can clearly observe the enhancement in 
the electronic-population change due to resonance effects. 
The very different picture that is obtained within the two 
different approximations emphasizes the importance of 
accounting for excitonic effects (also) in the strong field 
regime. 



V. SUMMARY 

We presented a novel approach to the ab-initio calcu- 
lation of optical properties in bulk materials and nano- 
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structures that uses a time-dependent extension of the 
BSE. The proposed approach combines the flexibihty of 
a real-time approach with the strength of MBPT in cap- 
turing electron-correlation. It allows to perform compu- 
tationally feasable simulations beyond the linear regime 
of e.g. second- and third-harmonic generation, four- 
wave mixing, Fourier spectroscopy or pump-probe ex- 
periments. Furthermore, being the approach based on 
the non-equilibrium Green's Function theory, it is possi- 
ble to include effects such as lifetimes, electron-electron 
scattering21 and electron-phonon coupling^ in a system- 
atic way. Finally, we have applied the TD-BSE to the 
case of /i-BN. First, we have calculated the optical ab- 
sorption and compared it with the results from equilib- 
rium ^i-MBPT validating our approach and numerical 
implementation. Then, we have shown the potentialities 
of the TD-BSE approach beyond the linear- regime by cal- 
culating the change in the electronic population due to 
the interaction with a strong quasi-monochromatic laser 
field. 
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Appendix A: An efficient method to update the 
COHSEX self— energy during the time evolution 

In this appendix we show how we store and update 
the 1]'^°'^*'°'' self-energy in a efficient manner. First of 
all we neglect the variation of the screened interaction 



W(r, r'; G<(t)) with respect to the G'^(r, r', t) by setting 
to zero the functional derivative dW/dG (see Sec. Ill Cp . 
Within this approximation the does not contribute 
to the time evolution, therefore only T,^°^ needs to be 
updated: 

S-''(r,r',t)=zI^(r,r') ^ ^„,k(r)'^:',k(r')G< „, JO- 

n,n k 

(Al) 

The KBE involves the matrix elements (to, k|S'"^^ |m', k) : 

s^r-'.kW= E p-,n(G')p:,'.„'(G)M/G,G'(q)G<„,(t), 

G,G',q ^'"^ k,q k-q 

(A2) 

where 

p™,„(G) = / ^:„,k(r)^„,k-q(r)e^(«+'i)"-. (A3) 

k,q J 

In order to rapidly update 'S^'^^ after a variation of 
G^(r, r', t), we store the matrix elements: 

M,n,m',n.n'^ E Pm,n (k, Q, G')p,*„/ (k, q, G)I^G,G' (q) > 

G,G' 

(A4) 

in such a way that T,'^^^,^, can be rewritten as 

^m!m',k(0 ~ E/ ^'^m,m' ,n,n' ' ^^_„' i^) ■ (^5) 
/ q.k k— fi 

q 

The M matrix can be very large, but its size can be 
reduced by noticing that: (i) the matrix M is Hermitian 
respect to the (to, to') indexes; (ii) the number of k and q 
points is reduced by applying the operation symmetries 
that arc left unaltered by the applied external field; (iii) 
for converging optical properties only the bands close to 
the gap are needed (see section ITV)) . As an additional 
numerical simplification we neglected all terms such that 
Mm,m\n,n'/ max{A/„ „ } < M^, whcrc Mc is a cutoff 

q,k q,k 

that, if chosen to be Mc — 5 ■ 10"'^ docs not appreciably 
affect the final results. In principle by using an auxiliary 
localized basis sel4S one can obtain a further reduction 
of the matrix dimensions, but in the present work we did 
not explore this strategy. 
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